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We report results from a systematic analytic strong-coupling expansion of the Bose-Hubbard model 
in one and two spatial dimensions. We obtain numerically exact results for the dispersion of single 
particle and single hole excitations in the Mott insulator. The boundary of the Mott phase can be 
determined with previously unattainable accuracy in one and two dimensions. In one dimension 
we observe the occurrence of reentrant behavior from the compressible to the insulating phase in a 
region close to the critical point which was conjectured in earlier work. Our calculation can be used 
as a benchmark for the development of new numerical techniques for strongly correlated systems. 

PACS numbers: 05.30. Jp, 05.70.Jk, 67.40.Db 
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Quantum phase transitions in strongly correlated sys- 
tems have attracted a lot of interest in recent years. In 
fcrmionic systems the Mott transition is complicated by 
the fact that in unfrustrated systems the antiferromag- 
netic transition and localization transition occur at the 
same point (see e.g. [[[)). For interacting Bose systems 
with spin zero the situation is much simpler and one can 
focus on the physics of the Mott transition. Strongly 
interacting bosonic systems are not only of academic in- 
terest. Physical realizations include Josephson junction 
arrays, granular and short-correlation-length supercon- 
ductors, flux-lattices in type-II superconductors and pos- 
sibly in the future ultracold atoms in a periodic potential 



To be specific we investigate the generic model for the 
Mott transition, the Hubbard model, for bosons (BH 
model) : 
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where the b\ and hi are bosonic creation and annihilation 
operators, hi = b\bi is the number of particles on site i, t 
the hopping matrix element, U > the on-site repulsion 
and \i the chemical potential. With short range interac- 
tions only the model has two phases at zero temperature: 
a superfiuid phase and a Mott phase. Much of the physics 
of the model was already understood qualitatively in an 
early paper by Fisher ct al. B and subsequent papers 
(see e.g. §f|,§). 

It is however interesting to obtain a quantitative un- 
derstanding of the model - for example to compare with 
experiments. To this end the BH model has been stud- 
ied numerically by Quantum Monte Carlo simulations 

in one and tw0 s P atial dimensions. 
Recently the one-dimensional case was also investigated 
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using the density-matrix renormalisation group (DMRG) 
fUSfl . This study found indications for an unexpected 
reentrant behavior from the superfiuid to the Mott- 
insulator as a function of the hopping amplitude t for 
certain values of the chemical potential. 

In this Letter we report for the first time a system- 
atic analytic strong coupling series to high order for the 
Bose-Hubbard model. Previous attempts which were re- 
stricted to rather low order fl7|| showed promising results 
but were not sufficient to investigate the asymptotic be- 
havior of the series. Recently M. Gelfand JL8| proposed 
a method for a linked cluster expansion with degenerate 
states. We have implemented the series expansion of the 
ground state and the first excited state as a linked cluster- 
expansion on a computer. The results show a spectacu- 
lar convergence of the Pade approximants for the phase 
diagram in one and two dimensions. The critical points 
can be determined to a previously unattainable accuracy 
(relative errors of w 10~ 3 ). In particular we are able to 
confirm convincingly that in one dimension there is reen- 
trant behavior of the Mott phase. The series calculation 
can be used as a benchmark for development of new nu- 
merical techniques for strongly correlated systems (e.g. 
DMRG). 

We start by writing down the ground state in the 
atomic limit (the hopping matrix element t — ► 0). In 
the atomic limit the number of bosons per site is fixed to 
an integer number say uq. Then the ground state of the 
Mott insulator with a fixed number no of particles per 
site is given by 
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with energy 
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Perturbation theory for the ground state energy Euott 
can be formulated as a linked cluster expansion, see e.g. 
p9| and the ground state energy can be obtained in the 
thermodynamic limit "relatively easily" . 

The Mott transition is obtained by studying charge 
excitations on top of the Mott phase. The charge excita- 
tions are gapped in the incompressible Mott phase and 
become gapless at the Mott transition. In the atomic 
limit charge excitations are created by adding or remov- 
ing a particle onto or from a particular site i: 
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Their energy relative to the ground state is given by 

fit/hole = ±(tf«o-/i) (6) 
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for particles and holes respectively showing that the 
charge excitations are degenerate. This degeneracy is 
lifted as soon as the hopping amplitude, t, is finite. In 
the atomic limit the energy of the excited states vanishes 
for a chemical potential Uc°^ = Uno and the system be- 
comes compressible. . 

A systematic strong coupling expansion of the energy 
of the charge excitations complicated due to the high 
degeneracy. The problem how to write down a linked 
cluster expansion for degenerate states was solved only 
recently by Gelfand |l8j . The idea is to construct pertur- 
batively an effective Hamiltonian Hf^ in the subspace of 

the degenerate states \rio; *)palt/hoio by a similarity trans- 
formation 



with #,,,(*) = S~}(t) 
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where Greek indices run over states in the full Hilbcrt 
space while Latin indices are restricted to the degener- 
ate manifold of single particle and single hole states (Q) 
and (^|) respectively. Then the linked cluster theorem ap- 
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plies to Hfj(t) — i?Mott(i). In the case of a homogeneous 



system Hfj depends only on the difference of indices 
i — j and is easily diagonalised by a Fourier transform. 
This way one can determine the full dispersion E(\&; t) of 
the charge excitations. In many ways the linked cluster 
expansion is similar to a exact diagonalization study of 
small systems - however in the linked cluster expansion 
it is possible to remove all finite size effects in each order 
and one obtains the full dispersion in the thermodynamic 
limit. 

For positive values of the hopping matrix element t the 
smallest (largest) eigenvalue in the particle (hole) sector 
is always located at a wavevector k = 0. The upper phase 
boundary of the Mott phase is thus determined by the 
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FIG. 1. Dispersion of the single-particle and single- hole 
excitations of the square lattice Bose-Hubbard model at 
t/U = 0.055. 



condition Bh i e (k = 0; t) — and the lower boundary by 
-E P art(k = 0;t) = 0. With increasing hopping t the dis- 
tance between the upper and lower boundary decreases 
until finally at some critical value, t c , the energy to re- 
move a particle and the energy to add a particle become 
degenerate and the Mott insulator vanishes altogether. 

We will first discuss the BH model, Eq. [j] on a two di- 
mensional lattice. We investigated both the square and 
triangular lattice and calculated the series for occupation 
numbers uq = 1 and no = 2 up to 12 th and 10 th order 
respectively. The dispersion of the particle and hole ex- 
citations for rig = 1 on the square lattice is shown in 
Fig.jl]. The different shape of the two curves reflects the 
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particle-hole asymmetry of the model Hamiltonian 
The series were found to converge very rapidly. Fig 
was obtained by summation of the 12 th order series. It 
turned out to be almost indistinguishable from the result 
of the 10-term series even for t/U = 0.055 which is not 
far from the critical endpoint t c of the Mott lobe. The 
particle and hole excitations both have a pronounced ex- 
tremum at wavevector k = and are separated by a gap 
A. For values of the chemical potential u in this range 
all single charge excitations are gapped and the system 
is insulating. 

The phase diagram shown in Fig. || is obtained by a 
Pade analysis of the series for the single particle gap, A. 
Scaling theory || predicts that in the neighborhood of 
the critical point (t c , u c ) the single particle gap A(t) as 
a function of the hopping matrix element, t, has the gen- 
eral form: A(t) = A(t)(t c — t) ZL/ where A(t) is a regular 
function of t and zv is the dynamical critical exponent. 
We use the following procedure to extrapolate the series 
p| . We calculate the logarithmic derivative of the series 
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FIG. 2. Phase diagram of the square lattice constructed 
using a Pade analysis of the series. The Mott phases are 
denoted by MI and the superfluid region by SF. The left 
curve is the lowest order Pade approximant (4 th order series) 
the right curve represents all the higher approximants. The 
inset shows a resolution of the region around the critical point. 
Note the scale! 



of the gap with respect to t which results in 

01og(A(*)) = Am () 

dt t-t e A{t) K ' 

The right hand side of Eq. || is well approximated by a 
Pade approximant. The pole of the Pade approximants 
for 91og(A(f))/9£ then determines the critical point t c 
and the residuum determines the dynamical critical ex- 
ponent zv. We then integrate the Pade approximants 
numerically to obtain the single particle gap A(t). With 
the exception of the lowest approximant all others ap- 
proximant turn out to be almost indistinguishable from 
each other indicating a rapid convergence. The results 
are shown in Fig. [2|. To observe any change at all in 
the higher approximants we have magnified the region 
around the critical point in the inset. The chemical po- 
tential is a regular function of the hopping matrix el- 
ement t. We used Pade analysis to check the scaling 
prediction and found for the critical point t c ~ 0.0599 
and the critical exponent v w 0.69. This has to be com- 
pared with the known value for the 3D xy- model pH , 
v = 0.6693 ± 0.0010, obtained by Borel summation of 
field theoretical results. The difference between the two 
results is of the order of a few per cent. Obviously the 
Pade analysis has a tendency to slightly overestimate the 
value of the critical point which in turn induces an error 
in the value of the critical exponent. 

It is also possible to extract the critical hopping ma- 
trix element, t, and the chemical potential at the critical 
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FIG. 3. 1/ order extrapolation of the critical point (i c ). 



point n c directly from the series. In each order k of the 
expansion the single particle gap A(i) vanishes at some 

(k) 

effective critical value t c with a corresponding effective 
fi{ k \ Plotting tc and vs. 1/k one finds again a 
rapid convergence as shown in Fig.(||). Extrapolation to 
k — » oo allows to determine accurately the critical point: 
t c = 0.05974 ± 0.00004 and fx c = 0.371 ± 0.001. 

We now turn to the one-dimensional case. From 
scaling theory || the critical behavior of the system 
is expected to be that of a Kosterlitz-Thouless transi- 
tion |2^] for which the gap closes according to A(t) cx 
A(t) exp (-W/V*KT - *) for |*kt-*| < 1, where A(t) is 
a regular function of t. The asymptotic form of the gap 
makes it difficult to approximate A(t) directly. Therefore 
we analyze the series for log(A(i)) 2 . The Pade analysis 
of the series yields spectacular agreement with the recent 
DMRG study of the phase diagram |l6| as is shown in 
Fig. U where we compare results from the series analy- 
sis with numerical data of QMC simulations by Batrouni 
and Scalettar [n| and DMRG data The agreement 
between the series and the DMRG data is excellent. Both 
calculations show that for a fixed chemical potential as a 
function of the hopping matrix element t the Mott phase 
is reentrant meaning that by increasing the kinetic en- 
ergy one returns to a localized state! The series analysis 
confirms the surprising behavior observed in the DMRG 
]l6| calculation. A simple intuitive way of understanding 
this phenomenon is the fact that Mott lobe is particle- 
hole asymmetric for the lattice problem. The Kosterlitz- 
Thouless behavior at the phase transition then implies 
a reentrant Mott phase. In higher dimensions this fea- 
ture can not exist because the transition has a power-law 
behavior. The accuracy of previous calculations was not 
sufficient to observe the reentrant behavior. The uncer- 
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FIG. 4. Comparison of the phase diagram obtained from 
series expansion (solid line), DMRG (solid circles) and QMC 
(solid squares). The Mott phase is denoted by MI and the 
superfluid phase by SF. 



taintics in the precise location of the Kosterlitz-Thouless 
transition are still comparatively large. We use a Pade 
analysis of In 2 A(i) oc (£kt — This quantity has a 

simple pole at the critical point which can be captured 
by rational function. This methods turned out to give 
excellent results. We estimate the point for Kosterlitz- 
Thouless transition to be located at tyr/U = 0.26 ± 0.01 
and hkt/U = 0.16 ±0.01. 

In real systems disorder plays an important role. With 
our method it is still possible to determine boundary of 
the Mott phase in that case by asking where the system 
becomes compressible. For a detailed discussion we refer 
the reader to Ref. [jl7|. The "Bose-glass" phase can not 
be studied with the techniques used above - since the 
groundstate has no gap. 

In conclusion, series expansion techniques were ap- 
plied to investigate the zero temperature properties of 
the Bose-Hubbard model in one and two dimensions. We 
determine the complete spectrum of single-particle and 
single-hole excitations in the Mott phase. The phase dia- 
gram in one and two dimensions is obtained quantitatively 
and the critical end points of the Mott insulator regions 
are determined. In two dimensions this is so far the only 
quantitative investigation of the complete phase diagram 
of this problem. In one dimensions the series shows al- 
most perfect agreement with a recent DMRG study and 
provides a conclusive confirmation for counterintuitive 
reentrance behavior from the compressible to the insu- 
lating phase near the Kosterlitz-Thouless point. 

We acknowledge useful and interesting discussions on 
this problem with M. P. Gelfand, T. Giamarchi, A. J. Mil- 
lis, A. v. Otterlo, R. R. P. Singh, G. Schon, H. Schulz. 
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